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ABSTRACT 

We present a model of magnetohydrodynamic (MHD) turbulence in the extended solar corona that 
contains the effects of collisionless dissipation and anisotropic particle heating. Recent observations 
have shown that preferential heating and acceleration of positive ions occurs in the first few solar radii 
of the high-speed solar wind. Measurements made by the Ultraviolet Coronagraph Spectrometer aboard 
SOHO have revived interest in the idea that ions are energized by the dissipation of ion cyclotron res- 
onant waves, but such high-frequency (i.e., small wavelength) fluctuations have not been observed. A 
turbulent cascade is one possible way of generating small-scale fluctuations from a pre-existing popula- 
tion of low-frequency MHD waves. We model this cascade as a combination of advection and diffusion 
in wavenumber space. The dominant spectral transfer occurs in the direction perpendicular to the back- 
ground magnetic field. As expected from earlier models, this leads to a highly anisotropic fluctuation 
spectrum with a rapidly decaying tail in the parallel wavenumber direction. The wave power that de- 
cays to high enough frequencies to become ion cyclotron resonant depends on the relative strengths of 
advection and diffusion in the cascade. For the most realistic values of these parameters, though, there 
is insufficient power to heat protons and heavy ions. The dominant oblique fluctuations (with dispersion 
properties of kinetic Alfven waves) undergo Landau damping, which implies strong parallel electron 
heating. We discuss the probable nonlinear evolution of the electron velocity distributions into parallel 
beams and discrete phase-space holes (similar to those seen in the terrestrial magnetosphere) which can 
possibly heat protons via stochastic interactions. 

Subject headings: MHD — plasmas — solar wind — Sun: corona — turbulence — waves 



1. Introduction 



In order to produce the hot (10 6 K) solar corona, a 
fraction of the mechanical energy in the Sun's internal 
convective motions must be converted into heat above 
the photosphere. Despite more than a half century of in- 
vestigation, though, the precise physical processes that 
lead to coronal heating and the subsequent acceleration 
of the solar wind are not known. At the coronal base, 
different combinations of mechanisms (e.g., magnetic 
reconnection, turbulence, wave dissipation, and plasma 
instabilities) are probably responsible for the varied ap- 



pearance of coronal holes, quiet regions, loops, and X- 
ray bright points (Priest et aL 2000). In the open mag- 
netic flux tubes that feed the solar wind, additional heat- 
ing at heliocentric distances greater than about 2 solar 
radii (Rq) is believed to be needed (e.g., Leer, Holzer, & 
FTa 1982; Parker 1991). This paper investigates the con- 
sequences of a promising mechanism for extended coro- 
nal heating: the kinetic dissipation of driven magnetohy- 
drodynamic (MHD) turbulence in the strong background 
magnetic field of the accelerating solar wind. 
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The necessity for gradual energy deposition at large 
distances from the coronal base comes from three sets 
of empirical constraints (see also Cranmer 2002). First, 
pressure-driven models of the high-speed component 
(v > 750 km s _1 ) of the solar wind cannot be made 
consistent with the relatively low inferred electron tem- 
peratures in coronal holes (T e < 1.5 x 10 6 K) with- 
out additional energy or momentum deposition. Second, 
spacecraft in the interplanetary medium have measured 
radial gradients in proton and electron temperatures that 
are substantially shallower than predicted from pure adi- 
abatic expansion (Phillips et al. 1995; Richardson et al. 
1995). Similar measurements of radial growth of the pro- 
ton magnetic moment between the orbits of Mercury and 
the Earth (Schwartz & Marsch 1983; Marsch 1991) point 
to specific forms of gradual energy addition. Third, the 
Ultraviolet Coronagraph Spectrometer (UVCS) aboard 
the Solar and Heliospheric Observatory (SOHO) mea- 
sured extremely high heavy ion temperatures, faster bulk 
ion outflow compared to protons, and strong anisotropics 
(with Tj_ 3> T||) in ion velocity distributions in the ex- 
tended corona (Kohl et al. 1997, 1998, 1999; Noci et al. 
1997; Cranmer et al. 1999b; Giordano et al. 2000). 

The list of possible physical processes responsible 
for extended coronal heating is limited severely by the 
fact that Coulomb collisions are extremely weak above 
2 to 3 Rq and that ions receive more energy than elec- 
trons (with Tj on » T p > T e ). Also, most suggested 
mechanisms involve the transfer of energy from propa- 
gating fluctuations — such as waves, shocks, or turbulent 
eddies — to the particles. This broad consensus has arisen 
because the ultimate source of energy must be solar in 
origin, and thus it must somehow be transmitted out to 
the distances where the heating occurs (see, e.g., Holl- 
weg 1978; Tu & Marsch 1995). The most common wave 
damping mechanisms proposed for the coronal base (re- 
sistivity, viscosity, and thermal conductivity) seem to be 
ruled out in the extended corona because of the relative 
unimportance of collisions. 

Wave-particle interactions are natural alternatives to 
collisional processes and have been studied in a solar 
wind context for several decades (Barnes 1968; Toichi 
1971; Abraham-Shrauner & Feldman 1977; Hollweg & 
Turner 1978; Marsch, Goertz, & Richter 1982; Isenberg 
& Hollweg 1983; Hollweg 1986; Tu 1987, 1988; Ax- 
ford & McKenzie 1992). The SOHO observations dis- 



cussed above have given rise to a resurgence of inter- 
est in the collisionless dissipation of ion cyclotron waves 
as a potentially important mechanism in the accelera- 
tion region of the high-speed wind (e.g., McKenzie et al. 
1995; Tu & Marsch 1997, 2001; Hollweg 1999b; Li et 
al. 1999; Galinsky & Shevchenko 2000; Cranmer 2000, 
2001, 2002; Hollweg & Isenberg 2002; Vocks & Marsch 
2002). In the extended corona, the cyclotron or Larmor 
frequencies of positive ions range from 10 to 10 4 Hz, 
but the oscillation frequencies observed on the solar sur- 
face (dominated by convection) are typically ^0.01 Hz. 
Thus, any wave generation mechanism seems to require 
bridging a gap of many orders of magnitude in frequency 
or wavenumber. 

Axford & McKenzie (1992) suggested that high- 
frequency Alfven waves could be generated during 
small-scale reconnection events in the rapidly evolving 
supergranular network. These waves would propagate 
up through the corona until they reach the cyclotron res- 
onance radii of various ions, then damp over a very short 
distance (see also Schwartz, Feldman, & Gary 1981). Tu 
& Marsch (1997) and Marsch & Tu (1997) presented so- 
lar wind models based on the idea that a launched power 
spectrum becomes "swept up" and eroded by the slow 
radial decrease of the cyclotron frequency. The gen- 
eral scenario of cyclotron sweeping has been called into 
question by Cranmer (2000), Hollweg (2000), and Lea- 
mon et al. (2000), but its relative importance in relation 
to other physical processes has not yet been fully estab- 
lished (see also Tu & Marsch 2001). 

Alternatively, there have been numerous "local" wave 
generation scenarios proposed that involve the transfer of 
energy from other sources into the ion cyclotron mode 
at a range of heights in the acceleration region of the 
solar wind. Into this class of models fall the processes 
of MHD turbulent cascade, kinetic plasma instabilities 
(driven by non-Maxwellian velocity distributions or spa- 
tial gradients), or wave mode conversion driven by re- 
flection or refraction (e.g., Hollweg 1986; Matthaeus et 
al. 1999; Kaghashvili & Esser 2000; Markovskii 2001; 
Moran 2002). Most of these local wave generation mod- 
els rely on the Sun launching a sufficiently intense spec- 
trum of low-frequency (i.e., with periods longer than ~5 
min) Alfven or fast-mode MHD waves that do not damp 
within a solar radius of the surface. At present, how- 
ever, there are few empirical constraints on the gener- 
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ation mechanisms, exact propagation modes, or power 
levels of ion cyclotron waves in the corona. 

In this paper we study the consequences of 
anisotropic MHD turbulent cascade as a possible gen- 
eration mechanism for ion cyclotron waves in the ex- 
tended corona. We also examine the turbulent cascade 
as a source of other forms of energy — such as parallel 
electron beams and phase-space holes — that could also 
lead to ion heating in the corona. We model the cascade 
as a diffusion process in three-dimensional wavenum- 
ber space and recover the basic form of the Goldreich & 
Sridhar (1995) spectral anisotropy. This solution for the 
turbulent power spectrum is coupled to a general Alfven 
wave dispersion relation in order to compute the rela- 
tive amounts of heating given to protons and electrons. 
The dispersive properties of the highly oblique kinetic 
Alfven wave (KAW; see Stefant 1970; Hasegawa 1976; 
Lysak & Lotko 1996; Hollweg 1999a; Stasiewicz et al. 
2000) turn out to be key drivers of the plasma condi- 
tions of the extended corona. The KAW has been stud- 
ied in a coronal and solar wind context by Dobrowolny & 
Torricelli-Ciamponi (1985), Song (1996), Leamon et al. 
(1998, 1999, 2000), Shukla et al. (1999), and Voitenko 
& Goossens (2002). 

The remainder of this paper is organized as follows. 
In § 2 we develop a semianalytical model of anisotropic 
MHD turbulence by extending the wavenumber diffusion 
picture of Zhou & Matthaeus (1990) to three dimensions. 
In § 3, proton and electron heating rates that are consis- 
tent with the derived power spectrum are computed from 
a Vlasov-Maxwell kinetic dispersion relation. For the 
realistic case that KAWs heat electrons in the direction 
parallel to the magnetic field, § 4 follows a speculative 
chain of reasoning that leads to the generation of elec- 
tron phase-space holes (EPHs) that can heat protons and 
heavy ions via Coulomb-like collisions. Conclusions, re- 
maining questions, and implications for future spectro- 
scopic observations are summarized in § 5. 



2. Turbulent Cascade in the Extended Corona 

In this section we derive the Fourier power spec- 
trum of turbulent Alfvenic fluctuations as a function 
of the wavevector k in the low plasma-beta extended 
corona. We assume that the turbulence is driven at 



large scales by waves that originate at the Sun, and 
that some process (such as reflection) generates inward- 
propagating waves to excite a cascade (Matthaeus et al. 
1999). The precise generation mechanisms of the low- 
frequency waves at the solar surface are unimportant in 
the context of this paper (see, e.g., Roberts 2000). The 
outward-propagating waves enter the corona and trans- 
form, in part, to a mixed population of outward and in- 
ward modes. Only when inward and outward wave pack- 
ets are allowed to "collide" can nonlinear couplings to 
higher wavenumber occur. The resulting cascade is a 
time-steady transfer of energy from large to small spa- 
tial scales that is believed to progress by successive mag- 
netic reconnection in the presence of the strong coronal 
"guide field" (e.g., Matthaeus & Lamkin 1986; Lazarian 
& Vishniac 1999). We assume that the turbulence be- 
comes fully developed on time scales that are short com- 
pared to the bulk solar wind outflow and the geometrical 
expansion of open flux tubes. This allows us to model 
the turbulence as spatially homogeneous in a small vol- 
ume element (with constant density and magnetic field 
strength) in the extended corona. 



2.1. Definitions of Fluctuation Quantities 

The energy density of fluctuations in a magnetized 
plasma is described most generally as a sum of electric, 
magnetic, and kinetic components. For linear (i.e., small 
amplitude) fluctuations, the total fluctuation energy den- 
sity SU is given by 

8tt 8tt 

M m 

where angle brackets denote averages over times much 
longer than the fluctuation time scales, SE is the per- 
turbed electric field, and <5B is the perturbed magnetic 
field (e.g., Bernstein 1958). The velocity distribution 
function of particle species s is assumed to be the sum 
of an undisturbed, stationary (zero-order) state /o and a 
small, linear (first-order) perturbation Sf. The masses 
and velocities of the particles are denoted m s and v, 
and the integration above is taken over all particle ve- 
locities for each species. For simplicity, we assume a 
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fully ionized plasma composed only of protons and elec- 
trons (see § 5 for a discussion of the impact of heavy 
ions). The quantity SU above is constructed as the sim- 
plest second-order and positive-definite quantity that is 
conserved exactly as a consequence of the fully nonlin- 
ear Vlasov equation (Davidson 1983). 

If the zero-order proton and electron velocity distribu- 
tions are assumed to be Maxwellian, the kinetic energy 
density term in eq. (1) can be expressed as the sum of 
fluctuations in bulk velocity and density (often referred 
to as "kinetic" and "thermal" energies, respectively). We 
assume that these two contributions to the total kinetic 
energy are the dominant ones, even when /o departs from 
a Maxwellian form. Thus, the total energy density is 
given by 



E 
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(2) 



where n s and T s are the zero-order number densities 
and temperatures, Sv and 5n are the first-order veloc- 
ity and number density perturbations, and k B is Boltz- 
mann's constant. For a nonrelativistic plasma, the elec- 
tric field term is usually negligible in comparison with 
the other terms. Plasmas with sufficiently strong mag- 
netic fields tend to be dominated by incompressible mag- 
netohydrodynamic (MHD) fluctuations (see, e.g., Gold- 
stein, Roberts, & Matthaeus 1995), which exhibit energy 
density equipartition between their magnetic and proton 
velocity fluctuations, 
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+ (Sv z ) p w (Sv 2 ) p , (3) 



where p is the total mass density, approximately equal 
to its proton contribution m p n p . At the base of the so- 
lar corona, the total velocity amplitude (SU / ' p) x l 2 is be- 
lieved to be of order 20 to 40 km s -1 (e.g., Mariska, 
Feldman, & Doschek 1978; Chae, Schiihle, & Lemaire 
1998). This population of waves is expected to be mostly 
propagating upwards from the Sun, with a radial am- 
plitude dependence close to that predicted by the con- 
servation of wave action (Hollweg 1973; Jacques 1977; 
Banerjee et al. 1998; Esser et al. 1999). 

In this paper, we need to evaluate several of the terms 
in the total energy density SU independently. How- 
ever, the turbulence as a whole is describable in terms 



of the total fluctuation power spectrum W, defined as 
the energy density per unit volume in three-dimensional 
wavenumber space and scaled as follows 

y = J d 3 kW(k h k ± ,t) . (4) 

MHD turbulence is dominated by a cascade of fluctua- 
tion energy in the direction perpendicular to the back- 
ground magnetic field, so let us also define the reduced 
power spectrum quantity Wj_, 

/+oo 
dk\\ W(k h k ± ,t) , (5) 
-oo 

which is proportional to the total energy density per unit 
In k±, integrated over the parallel wavenumber ku. The 
quantity W± has units of velocity squared, and is related 
to the total energy density via 
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(6) 



Note that ku can be both positive (for outward prop- 
agating waves) and negative (for inward propagat- 
ing waves), but the perpendicular wavenumber k± is 
positive-definite. 

The energy density spectra defined above are ap- 
propriate for following the flow of the total energy in 
wavenumber space, but they do not necessarily track the 
dominant turbulent motions on all spatial scales. One 
important example is the case of large-fc^ kinetic Alfven 
waves (KAW) that have wavelengths larger than the 
mean electron Larmor radius but smaller than the mean 
proton Larmor radius. These fluctuations are "felt" by 
all electrons in a similar manner as longer-wavelength 
MHD waves. However, some protons execute Larmor 
gyrations on a spatial scale larger than that of the fluctu- 
ations, and thus they become uncoupled from the waves. 
(Only the slowest protons — in the core of their velocity 
distribution — have small enough Larmor orbits to remain 
coupled to the electrons and waves.) Thus, even though 
the massive protons carry the bulk of the kinetic energy, 
it is the electron term (5v 2 ) e that most accurately de- 
scribes the smallest-scale fluctuations. Let us then define 
v 2 ^ as the power spectrum of electron kinetic energy fluc- 
tuations in the perpendicular direction, which is defined 
in a similar manner as eq. (6), 
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where only the perpendicular components of the elec- 
tron velocity fluctuation are considered. 1 In the low-fcj^ 
MHD limit, where protons and electrons oscillate with 
the same velocity, v\ is equal to W±. In general, we 
define 

v 2 ± (k ± ) = cj> e (k ± )W ± (k±) (8) 

where (f> e describes the departure from ideal MHD en- 
ergy equipartition, and is equal to 1 in the MHD limit. 
Formally, cf) e is the ratio of kinetic energy density in per- 
pendicular electron motions to the total energy density of 
the fluctuations at a specified wavenumber. The inverse 
quantity <p~ 1 is also approximately equal to the number 
fraction of protons that remain coupled to the fluctua- 
tions at the specified wavenumber. In § 3.1 we derive the 
exact behavior of (j) e as a function of k±, but for simplic- 
ity we present the following approximation, 



1 + k\R\ 



(9) 



where R p is the mean proton gyroradius, defined as 
the ratio of the proton most-probable speed w p = 
(2ksTp I 'nip) 1 1 1 to the proton cyclotron frequency fl p = 
eB/nipC. The above approximate form for <j> e was moti- 
vated by a simple estimate of the fraction of protons that 
remain coupled to fc^-scale fluctuations; i.e., 



1 



Jo 



d v± f p (v±) 



d 2 v± f p (v±) 



(10) 

For a Maxwellian velocity distribution with an effective 
"core" defined by vj_ < V m - ix = Q p /k±, the above ex- 
pression yields a fraction of 1 in the limit k±R p <C 1 and 
a fraction of (k±R p )~ 2 in the limit k±R p 3> 1. Eq. (9) 
is a simple function that smoothly bridges both limits. 



2.2. Dominant Two-Dimensional Cascade 

We describe the cascade of energy in MHD turbu- 
lence as a combination of advection and diffusion in 
wavenumber space. Chandrasekhar (1943) developed 
the general statistical theory of representing stochas- 
tic processes — such as Brownian "random walks" — as 



effectively diffusive. Leith (1967) proposed a gen- 
eral advection-diffusion equation for the evolution of 
fluctuation power in isotropic hydrodynamic turbulence. 
Leith's chief assumption was that the spectral cascade 
term must only redistribute energy in wavenumber space 
and should not alter its total magnitude. Thus, the spec- 
tral cascade term was written as the divergence of a flux- 
like quantity, where the primary contribution to the flux 
is from nearby regions of wavenumber space (e.g., Kol- 
mogorov 1941; Obukhov 1941). Pao (1965), Eichler 
(1979), Tu, Pu, & Wei (1984), and Tu (1988) modeled 
the spectral transfer as pure advection by deriving dimen- 
sionally consistent flux quantities that are scalar func- 
tions of the local wave power and wavenumber. Zhou 
& Matthaeus (1990) applied Leith's (1967) cascade phe- 
nomenology to MHD turbulence and derived diffusion 
coefficients consistent with both the Kolmogorov (1941) 
and Kraichnan (1965) energy transfer rates. 

It has been known for several decades that MHD tur- 
bulence in the presence of a background magnetic field 
develops a strong wavenumber anisotropy. Both nu- 
merical simulations and analytic descriptions such as 
"reduced MHD" (RMHD) have indicated that cascade 
from large to small spatial scales proceeds mainly in the 
two-dimensional plane perpendicular to the background 
field B (e.g., Strauss 1976; Montgomery 1982; She- 
balin, Matthaeus, & Montgomery 1983; Matthaeus et al. 
1996; Goldreich & Sridhar 1995, 1997; Maron & Gol- 
dreich 2001; Bhattacharjee & Ng 2001; Cho, Lazarian, 
& Vishniac 2002). We apply the Zhou & Matthaeus 
(1990) wavenumber diffusion idea to this largely two- 
dimensional transfer of energy from small k± to large 
k±. 

We model the interplay of wave injection (at small 
k±), cascade, and damping (at large k±) with an 
advection-diffusion equation similar to those used by 
Zhou & Matthaeus (1990), Miller, LaRosa, & Moore 
(1996), and Stawicki, Gary, & Li (2001). The k r 
integrated fluctuation power is assumed to obey the fol- 
lowing transport equation, 



dW±_ 

at 



d_ 

dx 



1 



-j3W± + 7 



dW± 
dx 



'For kinetic Alfven waves, the parallel electron velocity fluctuation begins to dominate when k± becomes of the same order as 
the inverse ion inertial length (see Hollweg 1999a). However, we restrict the definition of v± to the perpendicular fluctuations in 
order to be consistent with anisotropic MHD turbulence theory. 
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+ S(k ± ) - D(k ± )W± , (11) 

where x = In k±, t s is a characteristic spectral transfer 
time, and (3 and 7 are dimensionless constants. The 5* 
term describes the source of energy at low k±, which pre- 
sumably propagates up from the Sun in the form of low- 
frequency Alfven waves. The D term describes the dis- 
sipation of energy at high k± . The local spectral transfer 
time t s is assumed to be equal to the nonlinear time scale 
(k^v^)' 1 (Batchelor 1953; Goldreich & Sridhar 1995). 
We use the perpendicular electron spectrum v±(k±) to 
track the turbulent motions on all scales, which intro- 
duces a factor of 4> l J 2 into the advection-diffusion equa- 
tion when t s is written in terms of Wj_. 

Note that we did not include in eq. (11) any large- 
scale transport terms that depend on radial gradients of 
the coronal plasma parameters. This is in accord with 
the assumption discussed above that the cascade acts on 
a time scale that is short in comparison to representative 
flow times (due to wind advection, wave propagation, or 
geometrical expansion) over a coronal scale height. The 
interaction between large-scale "sweeping" terms in the 
corona and local turbulent cascade effects is discussed in 
detail by Leamon et al. (2000). 

The nonlinear cascade process modeled in eq. (1 1) re- 
flects some of the physics of wave dispersion and damp- 
ing, but not all. The spectral transfer time t s behaves 
differently in the ideal MHD regime (cj) e s=a 1) and in 
the KAW regime (<j> e oc k\). Moreover, the dissipa- 
tion term D in these regimes behaves as a power-law 
in k±, thus producing a phenomenology very similar 
to that of viscosity in fluid turbulence (e.g., Pao 1965; 
Leith 1967). However, we do not include high-frequency 
(proton cyclotron) dispersive effects in the cascade term, 
even though we do include them in the dissipation term. 2 
Such effects in the cascade term are expected to be neg- 
ligible because of the very limited extent of the "disper- 
sion range" in k\\ (see the closely-spaced solid contours 
in Figure 4a, below). Therefore, cyclotron dispersion 
cannot substantially affect the rate of energy transport 
from low to high fen. In future work we intend to ex- 
plore more self-consistent cascade phenomenologies for 
plasmas of arbitrary k\\ and k±. 



In order to specify the (3 and 7 constants as indepen- 
dent advection and diffusion strengths, we follow the sta- 
tistical analysis of van Ballegooijen (1986). In the "ran- 
dom walk" limit that individual turbulent displacements 
are small compared to the outer-scale correlation length, 
van Ballegooijen (1986) found that (3 « 7. We assume 
for simplicity that (3 and 7 are constants (i.e., indepen- 
dent of wavenumber). We also assume that (3 and 7 are 
of order unity because the primary physics of the cas- 
cade is encapsulated in the spectral transfer time t s . Our 
baseline model thus assigns (3 = 7 = 1, and other models 
vary either (3 or 7 while keeping the other equal to 1 . The 
most appropriate values of these constants to use in the 
extended solar corona should be determined ultimately 
from the analysis of numerical simulations or laboratory 
experiments. 

Figure 1 shows representative time-steady solutions 
to eq. (11) computed using a Crank-Nicholson numeri- 
cal finite differencing scheme. The Appendix contains 
further details about the solution method and the defini- 
tion of supplementary quantities. The low-wavenumber 
source term S(k±) was defined as a narrow Gaussian 
function that injects most of its energy at an outer- 
scale length L oul <~ fc~ t that we specified as approx- 
imately 3800 km. This value is consistent with the 
photospheric granulation scale (1000 km) expanded su- 
perradially in a polar coronal hole out to a distance 
of 2 R e (i.e., L out oc B Q l/2 ; see Hollweg 1986). 
The coronal hole magnetic field structure is assumed 
to follow the dipole/quadrupole/current-sheet model of 
Banaszkiewicz, Axford, & McKenzie (1998). Note 
also that a similar outer-scale length (~4200 km at 
the base) was derived — independently of granulation 
observations — by Chae et al. (1998), who balanced an 
empirical turbulent heating rate by local radiative losses 
in the transition region and low corona. We assume that 
the wave energy at the outer scale is generated at the solar 
surface, propagates outward into the corona, and under- 
goes some degree of transformation into a mixed popu- 
lation of outward and inward modes (by, e.g., reflection). 
The dissipation term D(k±) used in the models shown in 



2 One can avoid possible inconsistencies of this kind by assuming the wavelength ranges of cascade and dissipation are separated 
from one another. However, this allows only the computation of the total energy fluxes and not a detailed comparison of relative 
proton/electron and parallel/perpendicular heating rates; see also Leamon et al. (1999, 2000). 
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Figure 1 was assumed to be similar to that of standard 
kinematic viscosity; i.e., D oc k\. We plot the power 
spectra computed with three different normalizing con- 
stants for the dissipation term. 

The solutions for W±(k±) fall into several regimes 
of wavenumber. For k± <C fc ut> there is no signif- 
icant turbulent cascade, and the residual power decays 
like W± oc k^ 1 as a result of the ever-present diffusion. 
The power spectrum peaks around fc out at a maximum 
value that depends on the integrated total energy. In Fig- 
ure 1, we specify (SU/p) 1 ^ 2 = 50 km s _1 , which results 
in a maximum value of W 1 / 2 of about 10.2 km s _1 . The 
factor-of-5 difference between these two quantities is a 
simple consequence of the shape of W±(k±) around the 
peak at fc out . For k± > fc out , the power spectrum begins 
to decrease in an MHD power-law inertial range, with a 

—2/3 

Kolmogorov-Obukhov spectral index W± oc k ± (the 
standard one-dimensional power spectrum is given by 
W±/k± oc kj 5 ^ 3 ; see eq. [6] and the detailed deriva- 
tion in the Appendix). The spectral index in the iner- 
tial range does not depend on (3 or 7, but only on the 
wavenumber and power dependence of t s . As k± ap- 
proaches the inverse proton gyroradius, however, the fac- 
tor <f> e in t s becomes larger than 1, and the inertial range 
spectral index steepens from -2/3 to -4/3. We refer to 
these wavenumber regions below as the MHD and KAW 
inertial ranges, respectively (see also discussion of the 
"dispersion range" in interplanetary space by Stawicki et 
al. 2001). At high k±, where D becomes of the same 
order as the local cascade rate t~\ the power spectrum 
dissipates rapidly. 



2.3. Extension into Parallel Wavenumber 

The previous section described the dynamics of tur- 
bulence as a function of k± and ignored the behavior 
of the full power spectrum W as a function of fen. The 
strong predicted anisotropy of MHD turbulence justifies 
this approach to first order, but we are also concerned 
with the possible "leakage" of power to high values of 
fcii and thus to high frequencies. A substantial amount of 
power at the cyclotron frequencies of positive ions would 
lead to significant preferential ion heating. 

Let us model the spectral transport in fc|| with a dif- 
fusion term similar to that assumed above. It should be 



emphasized, however, that the actual parallel transport 
is probably nonlocal in wavenumber space (due to, e.g., 
intermittent high-order wave-wave interactions and non- 
linear steepening; see Medvedev 2000) and thus should 
not be describable by local diffusion. Our primary goal 
is to recover the phenomenology of the anisotropic "crit- 
ical balance" described by Goldreich & Sridhar (1995) 
and not to model the physics of nonlocal turbulent trans- 
port in detail. We thus propose the following advection- 
diffusion equation for W(k): 



dW 



_l d_ 

k± dk± 



{k ± v ± [-(3k\W- 



<9fc_L 



(k 2 ± W) 



+ ak_\_v±_k\\ 



'd 2 w' 

, dk 2 , 



+ S'-DW 
(12) 

which has been written in this way in order to recover 
eq. (11) when each term above is multiplied by k\ and 
integrated over fc|| . The dimensionless parameter a, pre- 
sumably of the same order as (3 and 7, determines the 
relative strength of the effective parallel cascade. The 
wavenumber scale is set by k\\, which is a typical paral- 
lel wavenumber determined from the condition of critical 
balance (i.e., strong coupling between turbulent mixing 
motions and Alfven-wave motions along the field). Let 
us define the dimensionless ratio 



y = 



k±V± 



= k \\l k \\ 



(13) 



where uj r is the real frequency of Alfven waves as a func- 
tion of fen and k±, and the Goldreich-Sridhar condition 
of critical balance is equivalent to the condition y = 1 . 
The source term <S"(k) in eq. (12) is related to the source 
term in eq. (1 1) as follows, 

/+00 
dk\\ S\k h k±) , (14) 
-00 

and we assume D(k) is dominated by the same kinetic 
sources of damping that were assumed in § 2.2. 

The advection-diffusion equation presented above 
does not depend on the sign of kn, and thus it produces 
diffusion for waves propagating in both directions along 
the magnetic field. In this paper we assume that W is 
symmetric about k» = (i.e., that the net cross helic- 
ity of the fluctuations is zero). This is consistent with the 
analysis of Goldreich & Sridhar (1995), but is manifestly 
not true for the lowest-A:^ (energy-containing) modes at 
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the outer scale in the corona. For these large-scale waves, 
there should be more power on the outward-propagating 
side (fc|| > 0) than on the inward-propagating side. How- 
ever, as the cascade proceeds to higher wavenumbers, 
the power in outward and inward propagating modes 
should become increasingly mixed (see, e.g., Roberts et 
al. 1992; Vasquez, Markovskii, & Hollweg 2002). The 
main reason we are modeling W(k) in this paper is to 
compute the kinetic dissipation and particle heating at 
high wavenumber, so the assumption that the outward 
and inward power is equal seems justifiable. 



In order to evaluate the quantity y defined above, 
we estimate the dispersion relation for kinetic Alfven 

1 /2 

waves (KAW) as oj r ~ k\\VA<f> e ■ This relation as- 
sumes frequencies much lower than the proton cyclotron 
frequency (see § 3.1 for a more complete discussion of 

1 12 

dispersion). The multiplicative factor of </> e represents 
the lower inertia in the KAW regime due to the decou- 
pling of protons from the wave motions. The coupled- 
proton fraction l/(j> e can be considered to be propor- 
tional to the effective "mass" to of a body undergoing 
simple harmonic motion; the frequency of oscillations 
for such a body scales as oj 2 <~ k/to, where n is the ef- 
fective "spring constant" determined here by the Alfven 

speed Va- Because both the numerator and denominator 

l li 

in eq. (13) contain factors of </> e , the quantity y can be 
expressed as 

hV A 

which is essentially the inverse of the Goldreich & Srid- 
har (1995) nonlinearity parameter (. We use this approx- 
imation to derive an analytic solution for W(k\i , k±), but 
in § 3.2 we replace the above relation with the complete 
definition for y (eq. [13]) and use a numerical solution 
for uo r . 

In the MHD inertial range, W± is proportional to 



k ± 2 ^ 3 , and thus y oc k\\k ± ^ ,J . Turbulent eddies that 
conform to the condition of critical balance (y = 1) ex- 



2/3 



2/3 



which im- 



hibit a wavenumber scaling of fc|| oc k 
plies they become increasingly elongated in the direction 
of the magnetic field as the cascade proceeds to higher 
wavenumber. This spectral anisotropy is discussed in 
more detail by, e.g., Goldreich & Sridhar (1995), Maron 
& Goldreich (2001), and Cho et al. (2002), and possi- 
ble observational evidence of such elongated inhomo- 
geneities was presented by Grail et al. (1997). In the 



KAW inertial range, the elongation should grow slightly 

— 1/3 

stronger, with y oc kiik ± 1 . However, dissipation is be- 
lieved to set in around k±R p <~ 1, before too much of 
the KAW inertial range can develop. 

To solve eq. (12) we make the important assumption 
that W is a separable function of two variables: k± and 
y. Following Goldreich & Sridhar (1995), we write the 
wave power spectrum as 



W(k ± ,y) 



V A W^ 



1/2 



(16) 



which was defined in this way to ensure that the dimen- 
sionless function g(y) is normalized to unity, 



f 

J —i 



dy g(y) = 1 



(17) 



and also that eqs. (4)-(6) are satisfied. The solution for 
W±(k±) is discussed above in § 2.2, and the Appendix 
shows how a known solution for W± allows eq. (12) to 
be approximated as an ordinary differential equation for 
g(y). The solution to this differential equation in the in- 
ertial range is 



g(y) oc 



7d - OV 



1 + 



a 



(18) 



with the normalization determined by eq. (17), and 
(/?/7) + C + l 



n 



2(1 - 



(19) 



(see Appendix for details). The dimensionless exponent 
£ is defined by 



2C 



k± dW± 
'W± dk ± 



(20) 



and in the MHD and KAW inertial ranges, £ = 1/3 and 
2/3, respectively. The above solution for g(y) superfi- 
cially resembles a generalized Lorentzian, or kappa dis- 
tribution (Olbert 1969), which is Gaussian for small ar- 
guments and approaches a power-law tail for large ar- 
guments. Cho et al. (2002) analyzed numerical turbu- 
lence simulations and found that g(y) may be fit reason- 
ably well by decaying exponential or Castaing functions. 
However, the simulations do not seem to have sufficient 
dynamic range in the high-fc|| modes to resolve a possible 
power-law tail as in eq. (18). 
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The relative amount of wave dissipation that heats 
protons and heavy ions depends sensitively on the power 
that diffuses to high frequencies; i.e., y ^> 1. In this 
limit, the above analytic solution in the inertial range 
behaves as g oc y~ 2n , which (after normalization) no 
longer depends on the parallel diffusion coefficient a. 
Thus, even when y » 1, the dominant diffusion still 
comes from the perpendicular j3 and 7 terms. The major 
contribution of the a term is to "bootstrap" the trans- 
port of power from low to high k\\; the perpendicular 
transport then takes over to determine the residual kn- 
dependence of the power spectrum. A key parameter in 
the y-dependence of the power-law tail is the ratio (5/ '7, 
which is the relative strength of perpendicular advection 
over that of diffusion. We extend our baseline assump- 
tion of the previous section and assume a = [3 = 7 = 1, 
but for cases when the ratio (3/-f is varied we keep the 
two diffusion coefficients equal to one another (a = 7). 



3. Wave Dispersion and Dissipation 

3.1. Vlasov-Maxwell Kinetic Theory 

Once the total fluctuation power spectrum W(k) is 
known, it becomes possible to compute the relevant par- 
ticle heating rates in the extended corona. Despite the 
fact that the strong, critically balanced turbulence de- 
scribed above is not very "wavelike" (i.e., a coherent 
wave survives for only about one period before nonlinear 
processes transfer its energy to smaller scales), we fol- 
low the standard procedure of using the linear wave dis- 
persion relation to compute frequencies and dissipation 
rates (see also Miller et al. 1996; Quataert 1998; Leamon 
et al. 1999). The amplitudes of the waves in the dissipa- 
tion range are extremely small (i.e., {SB 2 ) <C £?,,), so the 
assumption of linearity is probably adequate. In at least 
one other turbulent astrophysical plasma — the interstel- 
lar medium — a comprehensive study of various dissi- 
pation mechanisms (Spangler 1991) showed that linear 
damping processes are most likely to be stronger than 
nonlinear damping processes. 

We solve the linear dispersion relation for the real and 
imaginary parts of the frequency in the solar- wind frame 
(w = uj r + iu>i) assuming a known real wavevector k. The 



Vlasov-Maxwell dispersion relation, 

c 2 

-^[kx(kx E)] + e • E = , (21) 

is derived by linearizing and Fourier transforming 
Maxwell's equations (see, e.g., Stix 1992; Krauss- 
Varban, Omidi, & Quest 1994; Brambilla 1998). Above, 
E is the perturbed electric field, c is the speed of light, 
and e is the dielectric tensor constrained by the linearized 
Vlasov equation. The dispersion relation is essentially a 
3x3 matrix multiplying E, and we solve for the com- 
plex values of u> that make the determinant of this ma- 
trix zero. We use the numerical Newton-Raphson tech- 
nique to isolate individual solutions from a grid of start- 
ing guesses in u> r , Ui space. Besides the assumption of 
Maxwellian distributions (see below) and the use of nu- 
merical expansions for Bessel functions and the plasma 
dispersion function (Fried & Conte 1961; Poppe & Wi- 
jers 1990), the kinetic dispersion tensor is calculated ex- 
actly. The nine elements of e at each solution allow us 
to compute the electric field polarization vectors and the 
relative amounts of wave energy in electric, magnetic, ki- 
netic, and thermal perturbations for each wave mode (see 
eq. [2] and Appendix 1 of Krauss-Varban et al. 1994). 

We model the local plasma conditions at a heliocen- 
tric distance of 2 i? Q in the accelerating high-speed solar 
wind. The magnetic field is given by the polar coronal- 
hole flux tube model of Banaszkiewicz et al. (1998), with 
Sip = 7660 rad s _1 . The coronal-hole electron density n e 
at 2 Rq has been constrained to be between 2 x 10 5 and 
7 x 10 5 cm~ 3 , depending on the relative concentration of 
polar plumes (Cranmer et al. 1999b). We choose a value 
within this range that sets the Alfven speed Va to be ex- 
actly 2000 km s _1 . For the purpose of computing the 
wave dispersion, we assume isotropic Maxwellian distri- 
butions for the protons and electrons with temperatures 
constrained empirically as follows: T e = 8 x 10 5 K, 
T p = 2 x 10 6 K (Kohl et al. 1997, 1998; Esser et al. 1999; 
Doschek et al. 2001). The plasma beta, here denoted as 
f3 to distinguish it from the cascade advection strength (3, 
is defined as the ratio of the total plasma pressure to the 
magnetic pressure, 

^ = f*ZX*toT a « 0009 ' (22) 

s 

where the numerical value pertains to the coronal con- 
ditions described above. At all heights in the corona, 
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<C 1, but this quantity always remains larger than the 
electron-to-proton mass ratio m e /m p . 

In this paper we concentrate only on the solution 
branch corresponding to shear Alfven waves. This mode 
is expected to be dominant in the limit of nearly in- 
compressible MHD turbulence (Goldstein et al. 1995; 
Tu & Marsch 1995), and here we follow the evolution 
of Alfvenic fluctuations both to the high-frequency ion 
cyclotron limit and to the high-fc^ KAW limit. The 
Alfven mode is identified clearly at the lowest wavenum- 
bers in our computational grid (&||Va/^ p = 10~ 7 and 
k±R p = 10~ 6 ), and it is tracked incrementally as either 
fc|| or k± is increased slowly. We ignore the fast and slow 
magnetosonic wave modes — and thus possibly underes- 
timate the full heating rates in the extended corona — 
but these types of waves have long been thought to be 
damped substantially before they can reach the corona. 
We also ignore high-frequency plasma waves such as 
Langmuir and hybrid modes (which are expected to 
be much smaller in amplitude than the lower-frequency 
MHD waves) as well as other strongly damped solutions 
with no fluid counterparts (such as the Bernstein modes). 

Figures 2 and 3 show example solutions of the disper- 
sion relation for various ranges of wavenumber. Figure 
2a plots the real and imaginary parts of the frequency 
for a constant value of ku (in the low-frequency limit, 
uj r <C fi p ) as a function of k±. This figure thus illustrates 
the evolution of the KAW in the limit of large oblique- 
ness angles. The ratio \uji/u) r \ is negligibly small for 
small values of k±, but it grows to a value of order 0.2 in 
the KAW limit (kj_R p > 1). The damping for all plotted 
values of k± is dominated by the linear Landau reso- 
nance. Figure 2b plots the relative contributions to the 
total linearized energy density of the waves as defined 
in eq. (2). In the low-fc^ MHD limit, there is equiparti- 
tion between the magnetic field term and the proton ki- 
netic term (i.e., transverse velocity oscillations). When 
there is significant KAW dispersion, though, the elec- 
tron kinetic term becomes important as fewer protons re- 
main coupled to the small-scale fluctuations. "Thermal" 
density fluctuations also become non-negligible (Holl- 
weg 1999a). The protons and electrons exhibit identical 
relative density oscillations Sn/n, but the proton thermal 
energy density is larger because T p = 2.5T e . The quan- 
tity </> e (eq. [8]), which is defined as the ratio of kinetic 
energy in perpendicular electron fluctuations to the total 



energy density, is also plotted in Figure 2a. 

In the low-frequency, highly oblique KAW limit, 
the Alfven branch ceases to exist for perpendicular 
wavenumbers greater than k±R p ~ 14. This oc- 
curs when the parallel phase speed of the waves oj r /k\\ 
exceeds approximately two times the electron most- 
probable speed, and the number of electrons available 
to participate in coherent wave motions (in the tail of 
the distribution) becomes increasingly small. We believe 
this also explains the departure from the simple analytic 
KAW dispersion relation, 

uj r « k^V A ^\ + k\Rl , (23) 

where a temperature-dependent, order-unity factor mul- 
tiplying the k\ term has been omitted. The above ex- 
pression assumes full coupling between the waves and 
the entire velocity distribution. The cessation of solu- 
tions occurs at higher perpendicular wavenumbers for 
higher plasmas. 

For higher wavenumbers (k±R p > 14) in the low- 
plasma considered here, the kinetic extension of the 
fast-mode magnetosonic wave is the only surviving solu- 
tion with similar dispersion properties as the KAW. This 
mode is similarly longitudinal (i.e., with its electric field 
polarization vector nearly parallel to its wavevector), and 
it exhibits a similar distribution of energy density frac- 
tions as in the rightmost part of Figure 2b. Although we 
expect the turbulent cascade to deposit some power into 
this new mode once the KAW solutions disappear, we 
neglect this mode henceforth in order to determine how 
much of the turbulent energy is absorbed by the purely 
Alfvenic fluctuations. 

Figure 3 displays the dispersion properties of nearly 
parallel-propagating Alfven waves versus ku (with k±R p 
held constant at 10~ 3 ) as the wave frequency approaches 
the proton cyclotron frequency. The numerically com- 
puted frequency comes close to that of the cold-plasma 
approximation, 

Wr ~ k " VA ]f*~ir (24) 

(e.g., Dusenbery & Hollweg 1981). The damping rate Wj 
contains two contributions: one due to Landau damping 
(at low ku ) and one due to the proton cyclotron resonance 
(at high fc||). The waves in this regime are nearly incom- 
pressible, with the magnetic energy density varying from 
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equipartition with proton bulk motions in the MHD limit, 
to near dominance of the total fluctuation energy in the 
limit uu r — > Q p (see Figure 3b). 

We find that the Alfvenic proton cyclotron branch 
ceases to exist for parallel wavenumbers greater than 
knVA/^i p ~ 3.4. This "warm plasma" effect was pre- 
dicted to occur by Stix (1992) at about 



^||maxK4 



Va 



1/3 



(25) 



where wn p is the parallel proton most-probable speed, 
and this quantity is ^2.7 in the plasma conditions mod- 
eled here. The 20% difference between the predicted and 
actual cutoff values probably occurs because of the dif- 
ferent assumed proton and electron temperatures (elec- 
tron contributions are ignored in the above expression). 

We computed the real and imaginary parts of the 
Alfven-mode frequency u in a grid that extended from 
1(T 7 to 3.4 (in k\\V A /flp) and from 10" 6 to 14.2 (in 
k±R p ). Specifically, is the growth or damping rate 
(growth if positive; damping if negative) for the elec- 
tric field perturbation amplitude. In order to deter- 
mine the positive-definite damping rate for the domi- 
nant two-dimensional cascade of energy (proportional to 
the square of the amplitude), we set the quantity D in 
eq. (1 1) equal to -2o>j. We selected a locus of points in 
the wavenumber grid that traces out the critical balance 
condition of y = 1, with y defined preliminarily by the 
undamped spectrum W±. For most values of the wave 
power that we considered, the y = 1 curve occurs where 
the wavevectors are extremely oblique (k± ^> /c||), and 
thus they are in the KAW regime plotted in Figure 2. 
For these low-frequency fluctuations, the ratio LJi/u r de- 
pends only on k±, and we can write the frequency and 
the critically balanced damping rate along the y = 1 
curve as 



uv = k±v± = k ± 4> l J 2 w{' 



1/2 



D(k ± ) = ^2u; l = -^k ±( t>l /2 W 1 / 2 



(26) 



We found a parameterized fit for the k± dependence of 
the ratio uJi/ui r and combined it with eq. (9) to derive 
the following expression for the damping rate along the 



y = 1 curve: 



D(k±) w (5.95 x lO^cm" 1 ) 



1/2 



(1 + Q.lk\Rp) 1 / 2 ' 



(27) 

The numerical factor on the right-hand side applies 
specifically to the adopted coronal conditions at 2 i? Q . 
Also, this parameterization is valid only for the range 
of wavenumber where there are Alfvenic solutions (i.e., 
k±R p < 14). In the MHD limit, where k^R p < 1, 
D(kj_) is proportional to fc^ 3 , which is close to — but 
slightly steeper than — the commonly assumed viscous 
damping relation D oc k\. 



3.2. Quasi-linear Proton and Electron Heating 

The above calculation of the damping rate u>i does 
not reveal the relative amounts of energy absorbed by 
protons or electrons, nor does it describe how the parti- 
cle velocity distributions are driven away from isotropy. 
In this section we compute individual proton and elec- 
tron heating rates (in both the parallel and perpendicular 
directions) and bulk acceleration rates, based on assump- 
tions utilized by Quataert (1998), Quataert & Gruzinov 
(1999), and Marsch & Tu (2001). The individual heating 
and acceleration rates are computable in closed form in 
the quasi-linear limit (i.e., if <C 0J r )- We use this as- 
sumption to compute only the relative fractions of energy 
in the various rates, and normalize the total energy dis- 
sipation by the numerically computed (not quasi-linear) 
values of from the previous section. 

The approach used in this paper differs from that of 
Leamon et al. (1999), who isolated the electron and pro- 
ton contributions to the turbulent heating rate by compar- 
ing their baseline calculation with a trial calculation with 
an extremely low electron j3. The approach of Leamon 
et al. retained the full solution of the linear dispersion re- 
lation (even for \ui\ w u) r ), but at the expense of having 
to assume an unrealistically small electron temperature. 
Leamon et al. (1999) also evaluated the KAW dispersion 
relation in plasma conditions appropriate for the in situ 
solar wind (r > 20 Rq), whereas we are concerned with 
the corona. 

To retain the greatest generality, we assume a bi- 
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Maxwellian velocity distribution for each species s, 

n _ 



x exp 



w\ 



A. 



(28) 



where uu s is the bulk speed along the magnetic field and 
wu s and w± s are the anisotropic most-probable speeds. 
These quantities relate to the temperatures via 



wj s = 2k B T\\ s /m a 



w ±s 



2k B T ±s /m s . (29) 



Let us define the damping rate that is attributable to par- 
ticle species s as a positive quantity: 



,(w r -U|| s fc||)irf 



(0) 



(30) 



where uj ps = (4nqln s /m s y< 2 is the plasma frequency of 
species s and q s is the charge on particle s. The denom- 
inator of eq. (30) ensures that the sum of oji. s over all 
species must be equal to the magnitude of Wj. We define 
the generalized resonance function as 



(31) 



where lZ s (k, £) is the linear Landau/cyclotron resonance 
function given by eqs. (42)-(45) of Marsch & Tu (2001) 
and m is an integer. The quantity 1Z S depends on com- 
plicated functions of u>, k, and the electric field polariza- 
tion vector, and is too lengthy to reproduce here. Linear 
wave-particle resonances are included in this quantity as 
factors of the form 



exp 



(32) 



where I = corresponds to Landau damping and tran- 
sit time magnetic pumping (the latter believed to be im- 
portant only when j3 > 1), and t ^ corresponds to 
cyclotron or gyroresonant damping. In the models pre- 
sented below, we compute 1Z S fully for the Alfven modes 
discussed in the previous section, and we sum over 
I from —20 to +20. Test runs with larger ranges of I 
have ensured that no quantitative results depended on 
this choice of lower and upper summation bounds. 



Marsch & Tu (2001), following Melrose (1986) and 
others, wrote the particle heating rates as an integration 
over either a magnetic or an electric power spectrum. We 
have rewritten these rates in terms of W(k) and w, jS by 
noting that the total heating rate for species s (Q s ) should 
be balanced by the damping of the total fluctuation spec- 
trum, i.e., 



Qs = Q\\ s + Q±s = n s k B 



d_ (T\ 
dt 



•\\s 



+ T ±! 



(f k W(k) 2w t 



(33) 



The individual resonant acceleration and heating rates 
are thus given by an analogue of eqs. (27) and (33) of 
Marsch & Tu (2001): 



m s n s 



d_ 

dt 
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H 











p / k W(k) 




4(uj r 



(34) 



fey 

2n s K^/Kf^ 

For consistency with the assumptions in § 3.1, we have 
solved the above equations for isotropic Maxwellians 
(tU|| s = w± s ) and in the frame of the solar wind's bulk 
outflow (i.e., un s = 0). The resonant acceleration term 
presented above is given only for completeness. It has 
been shown that A\\ s is probably not an important con- 
tributor to solar wind acceleration, since the magnetic 
mirror force is typically stronger by a wide margin (e.g., 
Hollweg 1999b, Cranmer 2001). 

To evaluate eq. (34), the wave spectrum W{k\\ , k±) 
must be specified on the same grid of wavenumber points 
as in the dispersion calculation. The two-dimensional 
wave spectrum W±(k±) was first computed with no 
damping in order to specify the location of the y = 1 
curve, then it was recomputed with the damping rate 
D(k±) given in eq. (27). The change in W± did not 
significantly affect the position of the y = 1 curve in 
wavenumber space, so further iteration was not needed. 
We used the full numerical solutions for uj r and <j> e in 
the definition of y so that many of the simplifying as- 
sumptions in § 2.3 did not need to be applied. However, 
the analytic expression for g(y) in eqs. (18) and (19) was 
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used, with £ given by a simple undamped approximation 
that applies in the MHD and KAW limits, 



c 



1 / \+2k\R 2 p 
3 I 1 + k\R\ 



(35) 



The parallel diffusion strength a was always assumed to 
be equal to the perpendicular diffusion strength 7. The 
analytic expression for g(y) does not take account of cy- 
clotron damping at high values of y, which implies that 
the computed proton heating rates should be considered 
upper limits. 

Figure 4 displays various quantities in wavenumber 
space for the baseline model that assumes a = f3 = 7 = 1 
and (SU / p) 1 ^ 2 = 50 km s _1 . Contours of the proton and 
electron damping rates (dimensionless ratios |wj iS /w r |, 
for s = p, e) were computed in the same manner as the 
curves in Figures 2 and 3, separated into proton and elec- 
tron contributions using eq. (30), and plotted in Figure 
4a. In Figure 4a we also give indications of the angle 9 
between the background magnetic field and the wavevec- 
tor. Note that significant proton cyclotron damping oc- 
curs (for Maxwellian proton distributions) only when ku 
exceeds the inverse inertial length Q p /Va, with a rapid, 
Gaussian decline for smaller values of fcii. The onset 
of electron Landau damping, on the other hand, occurs 
much more gradually in the high-fcj^ KAW limit. 

Figure 4b shows contours of T4 / (k), computed as de- 
scribed above, and the critical balance (y = 1) curve. 
Although the decay of fluctuation energy as a function 
of fcii is steep, there is some residual energy that seems 
to make it to the proton cyclotron frequency. We should 
note, however, that for the case plotted in Figure 4b the 
solutions for k\\VA/Sl P > 10~ 2 are only approximate. 
The explanation for the power-law behavior of g(y) at 
large values of y depends on the existence of perpen- 
dicular wavenumber diffusion from high k± to low k±. 
As can be seen in Figure 4b, there is sufficient Alfvenic 
power available to be perpendicularly diffused only for 
fc|| Va/Q p < 10 -2 (i.e., to the right of the y = 1 curve), 
and the extension of the power law to larger values of fcy 
is an assumption that needs to be investigated further. We 
speculate that the inclusion of other (non-Alfven) modes 
at higher values of k± could make more power available 
to be diffused over a larger range of fc|| . 

The integration of eq. (34) provides the individual 
heating and acceleration rates. Below, we express these 



rates (i.e., Q|| s and Q± s ) dimensionlessly by dividing 
them by a fiducial empirical value Qemp = 10~ 8 erg 
cm~ 3 s _1 , which is representative of required heating 
rates at 2 R Q in various published models (e.g., Leer et 
al. 1982; Esser & Habbal 1995; Li et al. 1999; Dmitruk 
et al. 2002). For the baseline model illustrated in Figure 
4, the dominant dissipation channel is Landau damping, 
which is converted mainly into parallel electron heating. 
In this model, Q\\ e /Q e mp = 0.26, and Q± e is identically 
zero. There is a small amount of high-fc|| cyclotron ener- 
gization, which heats protons in the perpendicular direc- 
tion (Q± p /Q emp = 9. 1 x 10~ 5 ) and cools them slightly in 
the parallel direction (Q|| p /Qemp = —1-9 x 10~ 5 ). These 
rates, though, are far too small to provide the required 
energy to protons in the extended corona. 

We performed a number of rate calculations with 
different values for /?, 7, and SU. The total heating 
rates scale closely with the injected cascade rate £0 ~ 
(SU I p) 3 / 2 . Depending on the ratio (3/^f, though, as much 
as 30% of the input cascade energy can be unaccounted 
for in the total Alfven-wave damping. This is not surpris- 
ing because we cut off the wave spectrum abruptly when 
the Alfven branch solutions disappear and did not follow 
the flow of "mode-coupled" energy to higher wavenum- 
bers. Following the complete turbulent cascade of en- 
ergy on more than one dispersion branch is an important 
subject for future work. 

The following phenomenological scaling relations 
were produced as empirical fits to the numerical results 
found by varying (3, 7, and SU: 



Q\\e 
Qemp 

Q±p 
Qemp 



0.22 



(SU/p) 



1/2 



50 km s 



1/3 



(36) 



9.1 x 10~ 5 



(SU/p) 



1/2 



50 km s 



2n+l 



y^~ 2n ■ (37) 



The electron heating rate is relatively insensitive to the 
adopted values of f3 and 7, and its absolute value agrees 
with the empirically constrained total heating rate (i.e., 
Q|| e = Qem P ) if (SU/p) 1 / 2 w 120 km s" 1 . 

Note the sensitive dependence of the proton heating 
rate on the exponent n in the power-law tail of g(y). In 
the MHD inertial range, n is given simply by l+(3/3/47). 
The baseline case of (3 = 7 sets the normalization in the 
exponent above (i.e., for this case 2n = 3.5). The dimen- 
sionless factor y res , which we constrained by fitting to be 
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about 770, gives a rough indication of how far out in the 
tail the cyclotron interaction occurs. For proton heating, 
the dominant contribution to the integral in eq. (34) oc- 
curs when fc|| w k± (i.e., when k±R p w 0.2 at the proton 
cyclotron resonance) and at this wavenumber, 10 3 . 

Figure 5 plots the perpendicular proton heating rates 
from eq. (37) versus the ratio f3/j. Also plotted is a 
range of empirical proton heating rates from Li et al. 
(1999) that produce agreement with UVCS/SOHO H I 
Lya observations. It is clear that if (3 = 7, one would 
need unrealistically high values of the total wave power 
to agree with the empirical rates. However, the required 
levels of wave power are greatly reduced if 7 > (3. For 
example, Esser et al. (1999) derived an empirical range 
of 1 10 to 130 km s _1 for the total Alfven wave velocity 
amplitude at 2 Rq, using wave action conservation and 
spectroscopic limits on the base amplitude. To reach the 
empirical range of proton heating rates in Figure 5 simul- 
taneously with this constraint on (SU/p) 1 ^ 2 , one would 
need 7 to be greater than, or of the same order as, <~ 4/3. 
This does not seem to be an unreasonable condition on 
the turbulent cascade, but only simulations or laboratory 
experiments can provide firm limits on how the effective 
advection and diffusion strengths compare with one an- 
other. 

The use of SU as a total wave energy quantity glosses 
over the fact that, at the outer scale, there is probably 
considerably more power in outward-propagating waves 
(fc|| > 0) than in inward-propagating waves (fc|| < 0). 
Let us compare the heating rates computed above with 
the turbulence model of Matthaeus et al. (1999) and 
Dmitruk et al. (2001, 2002), who took the asymmetry in 
the outward and inward populations into account. Their 
net heating rate can be expressed as 



^out 

where Z_ and Z+ are wavenumber-integrated Elsasser 
amplitudes (with velocity units) of outward and inward 
propagating fluctuations, respectively, and L oul is a rep- 
resentative energy-containing length scale. Dmitruk et 
al. (2001, 2002) also approximated the reflected inward 
power Z+ as <~ L out \dVA/dr\. For the parameters at 2 
Rq discussed in § 3.1, the outer scale L out assumed in 
our cascade equation, and Q = Q smv , we solve for Z_ 
and Z + to be approximately 30 and 7 km s _1 , respec- 
tively. These quantities are amplitudes most comparable 



to the peak value of Wj_ at the outer scale, so in order to 
compare them with (SU / p) 1 ^ 2 we must multiply by the 
factor of 5 that arises because of the broadening of W± as 
a function of k± (see § 2.2). This implies a total outward 
wave power equivalent to (SU/p) 1 ^ 2 ~ 150 km s _1 , 
which is close to the empirical range reported by Esser et 
al. (1999). Our requirement above that (SU j p) x l 2 must 
be of order 120 km s~ 1 is also within this range. Thus, al- 
though we did not consider outward/inward asymmetries 
explicitly, the derived requirements on the total wave 
power seem robust and useful nonetheless. 

It is still uncertain whether the dissipation of Alfvenic 
turbulence in the extended corona will preferentially heat 
electrons or protons. If van Ballegooijen's (1986) sta- 
tistical limit of (3 = 7 holds for actual strong MHD or 
KAW turbulence, it seems clear that the electrons will 
be primarily heated along the magnetic field and there 
would not be sufficient wave power at the ion cyclotron 
frequencies to heat protons and heavy ions. This repre- 
sents a problem when confronted with spectroscopic ob- 
servations of extreme ion heating and anisotropy (with 
T± 3> T11) and at least mild preferential heating of pro- 
tons over electrons (e.g., Kohl et al. 1996, 1997, 1998; 
Esser et al. 1999; Cranmer et al. 1999b; Doschek et al. 
2001). We discuss alternate ways of energizing protons 
and heavy ions in the extended corona below in §§ 4-5. 



4. Nonlinear Evolution of Electron Beams 

Here we explore the possible consequences of the 
strong (but largely unobserved) parallel electron heat- 
ing that would be the result of Landau damping of ki- 
netic Alfven waves. The ideas discussed in this sec- 
tion are considerably more speculative than those dis- 
cussed above, but we are guided by detailed in situ ob- 
servations of roughly similar conditions in the Earth's 
ionosphere and magnetosphere. We gratefully acknowl- 
edge Matthaeus et al. (2003a, 2003b) for the initial sug- 
gestion that the following empirical constraints in near- 
Earth plasma environments could be applied to the solar 
corona and solar wind. 

A possible sequence of events could progress as fol- 
lows. The Landau damping of Alfven waves heats elec- 
trons in the parallel direction and has been shown, if al- 
lowed to develop far enough, to lead to non-Maxwellian 
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tails (Tanaka, Sato, & Hasegawa 1987, 1989; Leubner 
2001). These effectively beamed distributions may be 
unstable to the direct generation of ion cyclotron waves 
(Markovskii & Hollweg 2002) which could heat protons. 
However, we also note that electron velocity distribu- 
tions with monotonically decreasing tails may eventu- 
ally develop into non-monotonic bump-in-tail distribu- 
tions if the processes that accelerate electrons are in- 
termittent or otherwise spatially localized (e.g., Cairns 
1987). Even weak bump-in-tail electron beams have 
been shown to easily generate parallel Langmuir waves 
via various plasma instabilities (e.g., Thompson 1971; 
Melrose & Goldman 1987; Dum 1990; and references 
therein). Evolved Langmuir wave trains exhibit a peri- 
odic electric potential-well structure in which some of 
the beam electrons can become trapped. Adjacent po- 
tential wells can eventually merge with one another and 
gradually form isolated "holes" of saturated potential in 
electron phase space (Krasovsky, Matsumoto, & Omura 
1999; Omura et al. 2001; Mottez 2001). 

Electron phase-space holes (EPHs) are equilibrium 
electrostatic structures that have been predicted to exist 
for some time (Bernstein, Greene, & Kruskal 1957) and 
have been observed in the last decade in space plasma 
environments where electrons are accelerated along the 
background magnetic field (e.g., Matsumoto et al. 1994; 
Ergun et al. 1998, 1999; Bale et al. 1998, 2002). The tra- 
ditional formation mechanism for EPHs has long been 
considered to be the two-stream instability, but the ob- 
served EPHs have potential amplitudes that are too small 
to be formed from the strong counter-streaming elec- 
tron beams required to excite this instability. The above 
bump-in-tail/Langmuir excitation scenario is an alter- 
nate mechanism that has been shown to be able to form 
EPHs in numerical simulations. The stochastic nature 
of MHD turbulence (possibly dissipating in small re- 
connection sites) has also been shown to lead to similar 
kinds of intermittent nonlinearities (Dupree 1982; Am- 
brosiano et al. 1988; Biglari & Diamond 1989; Marsch, 
Tu, & Rosenbauer 1996; Tu, Marsch, & Rosenbauer 
1996; Drake et al. 2003). 

The reason we consider the effects of EPHs in the 
context of this paper is that in many places where phase- 
space holes are observed (e.g., Ergun et al. 1998), pro- 
tons exhibit preferential heating in the direction perpen- 
dicular to the magnetic field. Ergun et al. (1999) sug- 



gested that EPHs act as quasi-particles whose positively 
charged cores repel ions in a manner similar to ion-ion 
Coulomb collisions. Under certain conditions, a large 
number of these stochastic collisions can lead to ion 
heating. Because EPHs tend to flow along the magnetic 
field at velocities of order the electron most-probable 
speed (w|| e ~ 5000 km s _1 in the extended corona), the 
ions are roughly sitting still as the EPHs flow by. The net 
impulses the ions receive would thus be in the perpen- 
dicular direction. Below, we work out a quantitative es- 
timate for the perpendicular heating experienced by pro- 
tons in the extended corona due to EPH quasi-collisions. 



4.1. Phase Space Hole Properties 

Consider a homogeneous proton-electron plasma 
with the following localized disturbance in the electro- 
static potential: 



®(x, V, z) = $0 exp 



(39) 



where the z direction is that of the background mag- 
netic field and r 2 = x 2 + y 2 . The equipotential surfaces 
are oblate spheroids centered around the origin, and the 
above distribution is in agreement with several analyses 
of isolated EPH properties (e.g., Turikov 1984; Chen & 
Parks 2002). The parallel length scale zq has been ob- 
served in several space plasmas to be of the same order 
of magnitude as the electron Debye length, and we define 
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(40) 



where Ad 



5 /w pe \/2 is the electron Debye length 
and S is a dimensionless number that is observed to be 
typically 2 to 3. (Ergun et al. 1998). The transverse 
length scale tq may scale with the proton gyroradius in 
regions where ^l e /uo pe ^> 1, but Franz et al. (2000) esti- 
mated that 

,2 \ l ' 2 



ZQ 



1 + 



(41) 



in regions where fl e /u! pe <C 1 (which is more applicable 
to the solar corona). The above relation, which we apply 
in the model below, gives an aspect ratio r /zo of about 
3 to 10 in the extended corona. 
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The localized electric field resulting from the above 
potential can be determined analytically (E = — V$), 
and the charge density distribution is found by apply- 
ing Coulomb's law. The full analytic expression for the 
charge density p*(r, z) agrees with the plots in Figure 3b 
of Ergun et al. (1998), who modeled two extreme exam- 
ples of EPH geometry by effectively assuming r = z 
("sphere") and ro 3> zq ("plane"). There is a net positive 
charge density in the central core of the EPH which is 
surrounded by a net negative charge density. The posi- 
tive charge density at the center of the EPH (i.e., at r = 0, 
z = 0) is expressible in terms of the difference between 
local proton and electron number densities, 



Sn = n„ — n e = 



2ire 



2 1 

ri zi 



(42) 



which can be simplified further by neglecting the first 
term in the parentheses (because ro 3> zq). In this limit, 
we can express the ratio of Sn to the total electron den- 
sity 7l e as 

Sn _ 16 (e$ /m e ) 
n e S 2 is? 



(43) 



This allows us to convert between the central potential 
$o and the central charge separation Sn. 

We assume the EPH moves along the z direction with 
speed vu- This speed will probably be of the same or- 
der as the parallel electron thermal speed, so let us fol- 
low Turikov (1984) and define the Mach number M = 
va/w\\ e . Turikov (1984) found that one-dimensional 
EPHs become unstable and break up when M > 2, and 
Ergun et al. (1998) report values of M between ^0.2 
and 2 in the ionosphere. One cannot arbitrarily choose 
all three of the quantities 6, M, and Sn/n e . They are 
interrelated by the constraint that the electrons obey the 
Vlasov equation and can self-consistently affect the elec- 
tric potential via Poisson's equation. Turikov (1984) 
computed a relation that allows us to solve for 5n/n e 
as a function of S and M. For a Gaussian potential in the 
z direction (eq. [39]) and a weak charge separation (i.e., 
Sn/n e <C 1, which always seems to hold for reasonable 
values of the parameters), 



Sn 

n e 



-2M- 



(44) 



where Sq = 4ir 1 / 4 (2 In 4 — 1 ) 1 / 2 for the Gaussian poten- 
tial. Coincidentally, this value of So is extremely close to 



a value of 4.0. Chen & Parks (2001) noted that the above 
relation is formally only an upper limit on Sn/n e , but it 
is a useful scaling relation that we will assume holds for 
EPHs in the corona. Note from eqs. (43) and (44) that 
the central potential $o is proportional to S 6 , and thus is 
extremely sensitive to the parallel size of the EPH. 

One additional quantity that needs to be constrained 
is the spatial number density of EPHs — or equivalently, 
their filling factor — in the plasma. Only a detailed anal- 
ysis of EPH formation will yield constraints on how fre- 
quently they appear in space after significant merging 
has occurred, so here we use empirical limits from the 
in situ measurements. Defining the volume of an EPH as 
the equipotential surface that is 1 /e times its central po- 
tential, we relate their spatial number density n H to their 
fractional filling factor / H via 



n H 



h 



4wo/3 



(45) 



In a field of highly oblate EPH structures moving in 
the z direction, / H can be estimated by analyzing the 
mean time (or distance) between significant electric-field 
"kicks" felt by a slowly moving observer. For EPHs dis- 
tributed randomly in three-dimensional space, / H should 
be close to the ratio zq/Az, where the denominator is the 
mean distance between strong EPHs in the z direction. 
Ergun et al. (1998) found that there are often times when 
Az w vu/^lp, which we use in the rough calculations 
below. The specification for the filling factor also must 
take account of electron-electron Coulomb collisions, 
because if these collisions occur on a time scale faster 
than the time it takes to Landau-damp the MHD turbu- 
lence and build up electron beams, then EPHs should not 
form at all. Thus, let us modify the purely geometric fill- 
ing factor by an approximate collision term, with 
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(46) 



where Tdj ss is a turbulence dissipation time scale that we 
approximate as L out /Z_ (eq. [38]) and is about 200 s for 
(SU/p) 1 / 2 w 100 km s- 1 (see § 3.2). The standard elec- 
tron self-collision time is 



3/2 



Tcoll 



0.266 T t 
n e In A 



(47) 



(Spitzer 1962), where we assume In A = 21. When 
r coii *C T diss the EPH filling factor approaches zero, and 
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when Tcoii » rdi ss the filling factor approaches the value 
constrained by Ergun et al. (1998). For plasma condi- 
tions at 2 Rq, we can use canonical values of 5 = 3 and 
M = 1 to estimate / H to be ^10~ 4 , implying a spatial 
number density of 3 x 10" 



Tr z nn q\ $ 2 z % 



-10 cm -3^ 



4.2. Quasi-Collisions and Proton Heating 

A single electron phase-space hole moving along the 
z direction encounters a positive ion (with mass rrii and 
charge (/,) at a perpendicular impact parameter, or dis- 
tance of closest approach, denoted by r. We ignore 
Lorentz forces along the z direction because they are ex- 
pected to roughly cancel out between the approaching 
and receding halves of the EPH trajectory. The perpen- 
dicular momentum given to the ion over the entire trajec- 
tory is given by 



rrii Av±, 



/+C 
-ex 



dt qi E ± (t) 



(48) 



where the time t and the relative parallel displacement 
z between the EPH and the ion are related simply by 
z = vu t. The perpendicular electric field felt by the ion 
is that of the isolated EPH potential, 



E ± = (E x + Ey) 
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(49) 



We thus compute the perpendicular velocity impulse an- 
alytically to be 



2ir 1 / 2 q l r<S> Q z Q 

Avj_i = exp 

rmrQ v H 



(50) 



Following the standard development of Coulomb colli- 
sion rates (e.g., Spitzer 1962; Jackson 1975), we ex- 
tend the above calculation from one EPH to a "field" 
of EPHs distributed throughout space with number den- 
sity Tin- The number of ion/EPH encounters per unit 
time with impact parameters between r and r + dr is 
given by (2wr dr vhtlu)- The effective diffusion coef- 
ficient for perpendicular energization of the ion by the 
field of EPHs is then given by 



V 



dt 



f°° , 

/ dr2nr w H nn (Av±i) 
Jo 



to 2 vh 



(51) 



Because of the exponential dependence of E± as a func- 
tion of r, the above integral was well-defined and finite 
over all values of the impact parameter. This is a rel- 
ative simplification over the ion-ion Coulomb collision 
problem, where the above integral formally diverges and 
must be truncated using a maximum impact parameter. 

In the limit of a field of very low-energy ions, T>±i is 
essentially the heating rate per particle, per unit mass. As 
the ions heat up, however, the amount of energy they ex- 
tract from the EPHs cannot be arbitrarily large. Ideally, 
we should solve coupled energy conservation equations 
for both ions and EPHs, where each equation would con- 
tain energy exchange terms that would eventually lead to 
equipartition. It is not clear, though, how much of the 
EPH energy is actually "disposable" (i.e., easily given 
away to ions in collision-like interactions). 

In the frame at rest with respect to a phase-space hole, 
its electrostatic potential energy is given by the following 
volume integration, 
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(52) 



which is considerably larger than the kinetic energies of 
individual protons or electrons in the solar corona. In 
the frame of an ion, the EPH is moving by at a relatively 
high speed (though we assume v H <C c), and an induced 
perpendicular magnetic field is felt by the ion with a po- 
tential energy given by 
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which is several orders of magnitude smaller than the 
total electrostatic energy of the EPH. The induced mag- 
netic energy is proportional to and seems to be anal- 
ogous to a kind of kinetic energy. Our conjecture is that 
represents an order-of-magnitude estimate for the en- 
ergy available to be given away by EPHs in collisions 
with ions. We express this energy quantity as an equiv- 
alent temperature (Ub = ^b^b) and insert the adopted 
conditions at 2 Rq (i.e., A_d = 14 cm and i0|| e = 5700 
km s~') to obtain 



T B w (4.3 x 10 I2 K) ( j ) MV 4M2 (54) 
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(see, e.g., eqs. [43] and [44]). For increasing M, the to- 
tal magnetic energy in the EPH begins to decrease above 
about M rts 0.2 because very fast EPHs contain fewer 
trapped electrons and thus cannot support as high a po- 
tential as a slower EPH. 

In the extended corona, the value of Tb is huge in 
comparison to individual particle temperatures. Further, 
this large value of T B is smaller by at least a factor of 
v\j<? than the temperature quantity equivalent to the 
total electrostatic energy of a phase-space hole. EPHs 
seem to embody a substantial "pool" of energy that is 
continually replenished by KAW turbulence and Landau 
damping, and this energy does not seem to be diminished 
significantly by collisions with ions. Thus, we do not 
consider the approach to energy equipartition any fur- 
ther and assume that eq. (51) gives the effective net ion 
heating rate per unit mass. 

For the specific case of protons, we solve a steady- 
state solar wind internal energy conservation equation 
with the EPH diffusion coefficient as the sole source of 
extended heating. Our approach ignores the impact of 
heat conduction and energy sharing via particle-particle 
collisions, but these effects are expected to be relatively 
weak in the extended corona (e.g., Olsen & Leer 1996; 
Li 1999). We thus assume that the proton perpendicular 
temperature T± p evolves with radius according to 
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(55) 



where u\\ p (r) is the time-steady outflow speed profile and 
A(r) is the cross-sectional area of a solar wind flow tube. 
The second term in parentheses above is responsible for 
adiabatic proton cooling due to the geometric expansion 
of the plasma. Empirically determined values for the 
radial dependences of un p , A, n e , Bo (the background 
magnetic field strength), and T e (the electron tempera- 
ture; needed to compute A_d) are elaborated by Cranmer, 
Field, & Kohl (1999a). We integrate eq. (55) from an in- 
ner boundary at the base of the corona to an outer bound- 
ary in the extended corona. The remaining free param- 
eters are the inner boundary value of T± p and the EPH 
parameters 5 and M, which we assume to be constant in 
the extended corona. We assume the EPH filling factor 
/ H is given by eq. (46). 

Figure 6a shows the effective heating rate Q^ p = 
pT>j_ p /2, computed using eqs. (43^47) and eq. (51), for 



various values of S and a fixed value of M = 1. It 
also plots Q for the Dmitruk et al. (2001, 2002) turbu- 
lence dissipation rate (eq. [38]), with the radial depen- 
dence pZ 2 _ oc Va(u\\ p + Va)~ 2 A~ 1 given by wave ac- 
tion conservation (see also Zank, Matthaeus, & Smith 
1996). We interpret the order-of-magnitude similarity of 
the two sets of rate quantities as a promising plausibility 
argument for EPH proton heating. Despite this similar- 
ity, though, the two sets of rates are essentially "apples 
and oranges," since the Dmitruk et al. values come solely 
from the turbulent cascade rate and the EPH values do 
not contain any information about the cascade. The EPH 
curves in Figure 6a are not energetically consistent with 
the turbulence because we use the empirical filling fac- 
tor /h discussed above. A more self-consistent model of 
EPH heating would follow the turbulent energy budget 
from the dissipation by Landau damping to the growth 
of electron beams (which would embody only a fraction 
of the total dissipated energy) to the eventual intermittent 
formation of EPHs. 

For protons that start out with T± p = 10 6 K at r = 
1 .25 Rq (assuming some collisional contact with elec- 
trons), Figure 6b plots the computed radial dependence 
of T± p for the same range of S values as in Figure 6a. 
(The temperature curves are similar to those computed 
for a single value of 5 and a range of M values.) Sub- 
stantial proton heating is evident, but the results are very 
sensitive to the values of the EPH parameters. It would 
probably be more realistic to assume distributions of 5 
and M (see, e.g., Figure 4a of Ergun et al. 1998), but this 
introduces free parameters as well. Note, though, that the 
mean value of 5 reported by Ergun et al. (1998) is <~3.6, 
which is closest to the value that would give T± p (r) most 
similar to the UVCS H I empirical model results (Kohl 
et al. 1998). 

The proton heating models shown in Figure 6 are sug- 
gestive but nowhere near conclusive. A major source 
of uncertainty is the application of the ionospheric EPH 
spacing relation Az = vn/Ct p to the calculation of the 
EPH filling factor in the corona. This could be differ- 
ent by orders of magnitude, which would significantly 
change the above results. Also, our assumptions that S 
and M are constant with radius and that eq. (44) gives 
the exact central EPH potential should be examined crit- 
ically. More work should be devoted toward investigat- 
ing the growth and development of EPHs in the corona, 
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as well as their interactions with protons, heavy ions, and 
the "untrapped" electrons. 



5. Discussion and Conclusions 

In this paper we suggested two alternative "channels" 
by which protons and other positive ions in the extended 
solar corona can be heated perpendicularly to the back- 
ground magnetic field (see Figure 7). First, the parallel 
decay of anisotropic MHD turbulence to ion cyclotron 
frequencies may occur if the cascade is properly mod- 
eled by advection and diffusion strengths (3 and 7 with 
the property 7 > 4(3. This result could lend some sup- 
port to the many other proposed models that describe the 
damping of parallel-propagating ion cyclotron waves in a 
coronal context. (Note, however, that many more issues, 
such as implementing a more self-consistent cascade and 
damping formalism at large values of fc||, still needs to 
be resolved). Second, we proposed that, if 7 < 4(3, the 
dominant parallel electron heating associated with KAW 
Landau damping could lead to the nonlinear generation 
of electron phase-space holes. These EPHs have been 
shown to undergo collision-like interactions with pos- 
itive ions and heat them mainly perpendicularly. The 
EPH idea is extremely speculative because as yet we 
have no quantitative constraint on their spatial filling fac- 
tor in the corona. Only if they are generated in sufficient 
numbers (i.e., only if the conversion from KAWs to elec- 
tron beams to EPHs is efficient enough) can they provide 
substantial proton heating. 

There are several other related mechanisms by which 
protons and heavy ions could be heated in the extended 
corona. If there is substantial power in obliquely prop- 
agating fast-mode magnetosonic waves, their collision- 
less damping may contribute to proton heating (e.g., Li 
& Habbal 2001; Hollweg & Markovskii 2002). Oblique 
Alfven and fast-mode waves can steepen into shocks un- 
der certain conditions, and numerical simulations that 
employ the derivative nonlinear Schroedinger equation 
(DNLS) have produced a rich variety of steepening phe- 
nomena that produce power at high-frequency harmon- 
ics of an input spectrum (e.g., Spangler 1997). Certain 
types of collisionless shocks may also accelerate ions in 
the direction perpendicular to the magnetic field (Lee & 
Wu 2000). The basal generation and "sweeping" of high- 



frequency waves (e.g., Axford & McKenzie 1992; Tu & 
Marsch 1997) must also be developed further to deter- 
mine to what degree these waves can survive in the ex- 
tended corona (see also Cranmer 2001). Because of the 
multiplicity of time and spatial scales involved in study- 
ing many of these processes, the dominant physics may 
become evident only when a large number of competing 
mechanisms are included together in a solar wind model 
and allowed to evolve self-consistently. 

The proton heating that arises from ion cyclotron 
wave dissipation depends sensitively on the shape of 
the proton velocity distribution function. We assumed 
Maxwellian distributions in the above analysis, and the 
conclusions are not affected strongly if there is moderate 
bi-Maxwellian anisotropy. If the protons have suprather- 
mal high-energy tails (e.g., Scudder 1992), though, the 
resonance function 1Z S would have a less steep frequency 
dependence compared to eq. (32). In this case, waves 
with lower ku would be included in the thermally broad- 
ened ion cyclotron resonance and thus lead to a a lower 
effective value of y res in eq. (37). For strong enough 
suprathermal tails, it is possible that even the baseline 
cascade model (with f3 = 7) could contain enough 
parallel decay of wave energy to heat protons substan- 
tially. It should be noted, however, that Isenberg (2003) 
concluded from a completely collisionless description 
of "kinetic shell" velocity distributions that the damp- 
ing of parallel-propagating ion cyclotron waves is not 
a likely source of proton heating even if there is suffi- 
cient wave power in the extended corona. Indeed, damp- 
ing rates for the non-Maxwellian distributions predicted 
in quasi-linear models of collisionless resonant heating 
can be substantially different than those computed for 
Maxwellians (e.g., Isenberg 2001; Cranmer 2001; Tu & 
Marsch 2002). Full kinetic models that follow the pro- 
tons and electrons from the collisional low corona to the 
collisionless extended corona are needed to resolve these 
issues. 

Some mention should be made of the effects of heavy 
ions on the above analysis. It has been clear for sev- 
eral decades that the 5% to 10% abundance of He 2+ in 
the corona is responsible for a strong resonance in the 
Alfven wave dispersion relation, but it is unclear whether 
the other thousands of minor ion species are numerous 
enough to affect the dispersion relation (Isenberg 1984; 
Gomberoff, Gratton, & Gnavi 1996). Mode coupling be- 



-20- 



tween the different branches must be modeled explicitly 
in order to determine how much wave power survives at 
the proton cyclotron resonance. Also, it should be made 
clear that the EPH quasi-collision mechanism proposed 
above is not an efficient source of preferential heating for 
heavy ions, even if it is able to supply sufficient heat to 
the protons. The effective ion heating rate is proportional 
to rrii'D±i, which scales with charge and mass as qj/rrii 
(eq. [51]). The 5+ ions measured by UVCS/SOHO in 
the extended corona have a perpendicular kinetic temper- 
ature about two orders of magnitude larger than that of 
protons. The EPH mechanism, though, would only pro- 
vide an oxygen temperature larger than the proton tem- 
perature by a factor of 5 2 / 16 <~ 1 .6. The damping of ion 
cyclotron waves still seems to be the most likely source 
of preferential heating and acceleration for heavy ions. 

Improvements in remote-sensing measurements of 
the extended solar corona are needed to make signif- 
icant further progress in identifying and characterizing 
the most important physical processes. Next-generation 
spectroscopic instruments are being designed with the 
capability to sample the velocity distributions of dozens 
of ions in the acceleration region of the high-speed wind, 



as opposed to just 2-3 ions with UVCS. In addition to be- 
ing sensitive to many more emission lines, these instru- 
ments could also detect subtle departures from Gaussian 
line shapes that constrain the presence of specific non- 
Maxwellian distributions (e.g., Cranmer 1998, 2001). 
Seldom-used diagnostics such as the measurement of the 
Thomson scattered H I Lya profile (which probes the 
line-of-sight electron velocity distribution; see Withbroe 
et al. 1982) can put firm constraints on both the "core" 
electron temperature and the existence of beam-like or 
power-law wings. Improvements in radio sounding ob- 
servations are also making it possible to extract a great 
deal of information about MHD turbulence in the solar 
wind (e.g., Spangler 2002). 
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A. The Advection-Diffusion Cascade Model 



This Appendix provides details concerning the solution of eqs. (11) and (12). The two-dimensional advection- 
diffusion cascade model (eq. [11]) can be simplified by defining the following supplementary variables: 



W(k±) 



3/2 ,-3/3/2-7 



W^k ± 



We thus can define the perpendicular cascade rate e as 



-F(k ± ) 



dW 



which in the MHD regime (<j) e = 1) is proportional to k±v 3 ± . With the above definitions, eq. (1 1) simplifies to 

dW\ de 

= -k±— + S(k±) - D(k±)W± . 



(Al) 
(A2) 

(A3) 
(A4) 



dt ^ dkj_ 

In the inertial range, where both S and D are negligibly small, it is easy to see that the time-steady solution exhibits 
a constant value of s, which we will denote e . This cascade rate is supplied by the source function S(kj_), which we 
estimate to have the following Gaussian shape: 

2" 

/ X X all t } 

■ exp 



S(x) 



7T 1/2 0-out 



Cout 



(A5) 
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where x = lnk± and x out = In fc out . The dimensionless width of the Gaussian is given by er out = 1 in all models 
presented in this paper. With this explicit form for S, the wavenumber dependence of the time-steady cascade rate (at 
wavenumbers where D w 0) can be found by integrating eq. (A4), with 



six) = - 



1 +erf 



^out 



(A6) 



This provides a smooth transition from zero cascade (when k± <C fc out ) to cascade at the specified rate £o (when 
k± >• fc out ). Thus, the solution for W in the inertial range can be computed straightforwardly via 

~ f +ca e(k\ ) dk\ 

J k± Fik'J 

When k± <C k oat , the cascade rate approaches zero and W is constant, implying an energy spectrum W±_ oc k^ 1 . In 
the MHD inertial range, e is equal to a constant value of eq and eq. (A7) can be integrated analytically to yield that 

3/2 I 2/3 

Wjj is proportional to k ± . This results in the Kolmogorov-Obukhov spectrum W± oc k ± .In the KAW inertial 
range, <p e oc k\ and eq. (A7) integrates to obtain W± oc k ± 4 ^ 3 . 

For the full range of k± values below the dissipation range, the integration of eq. (A7) was performed numerically 
to obtain the undamped curve in Figure 1 . For 4> e given by eq. (9), an approximate "bridging solution" for the inertial 
range is 

3/T 



3eo ^ 2/3 r 



1+2^ I \/l + k]R 2 n + k^R, 



-2/3 

(A8) 



, 2 1 

which is valid in the MHD and KAW inertial-range limits and differs by no more than 6% from the exact numerical 
solution when k±R p w 1 . The above solution was used as the initial condition in the Crank-Nicholson finite difference 
code (e.g., Press et al. 1992) that included the effects of dissipation. All runs of the code reported in this paper utilized 
an evenly spaced grid in In k±, spanning six orders of magnitude from k±R p = 10~ 3 to 10 3 . The portions of the wave 
spectrum below the minimum value of k± were assumed to remain at the values given by the numerical integration 
of eq. (A7). The implicit nature of the finite-difference scheme allowed us to use a natural time step of 0.1 times 
the minimum value of t s in the wavenumber grid. The system was evolved toward a steady state, and any numerical 
diffusivity was identified by simultaneously evolving the same spectrum both with (D / 0) and without (D = 0) 
damping. The latter case isolated long-time-scale numerical effects that were negligible in magnitude, but nevertheless 
divided out of the final damped spectrum. 

An interesting feature of eqs. (11) and (A4) is that the behavior of the inertial range is largely independent of 
the values of (3 and 7. Indeed, even in the limit of pure advection (i.e., 7 = 0), one can derive analytic solutions for 
the inertial range and the dissipation range. (Some of the supplementary quantities defined in this Appendix become 
undefined for 7 = 0, but the original advection-diffusion equation remains well-defined.) For 7 = 0, the rapid cutoff in 
the dissipation range is an explicit exponential decline (see, e.g., Townsend 1951; Pao 1965). 

The remainder of this Appendix is devoted to the derivation of eq. (18), the analytic solution for g(y). Beginning 
with the time-steady advection-diffusion equation (eq. [12]) in the inertial range (S = D = 0), we substitute in the 
ansatz functional form for W(k±, y) given in eq. (16). Then, to evaluate the partial derivatives with respect to k±, 
we assume that W± conforms to a local power-law form proportional to k^ (as defined in eq. [20]). Even if £ 
is allowed to vary as a "slow" function of k±, the dominant variation in W± should remain the power-law decline 
with the local value of £. This assumption allows ( to be considered constant as a function of y. After considerable 
but straightforward algebra, eq. (12) becomes the following ordinary differential equation with y as the independent 
variable: 

(2C - V) Hy) + (1-0V h'(y) + a g"(y) = . (A9) 
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The supplementary quantity h is defined as 

h(y) = [/? + 7(i + 0]<7(y) + 7(i-C)y<?'(y) , (Aio) 

and primes denote differentiation with respect to y. The dimensionless exponent i] is defined as 

^ 2<j) e dk± 

(i.e., t) = in the MHD inertial range and i] = 1 in the KAW inertial range). Serendipitously, the ratio 

2C-V 

i-C 

is equal to 1 in both the MHD and KAW inertial ranges. From numerical solutions of W± and <f> e , the above ratio is 
never less than 1 or larger than ~1.3 in the narrow transition between the MHD and KAW inertial ranges. Thus, it 
seems reasonably safe to assume that this ratio is unity for all relevant wavenumbers, and we simplify eq. (A9) as 

(1 - [Hy) + yh'(y)] + ag"(y) = . (A12) 

This equation is integrated once, assuming g'(0) = 0, to obtain the first-order differential equation 

(l-Qyh(y) + ag'(y) = . (A13) 

Inserting the definition of h(y), we obtain the relatively simple equation 

Idg -(l-C)[/3 + 7(l+C)]2/ 



gdy a + 7(l-0 2 2/ 2 



(A14) 



which can be integrated analytically to obtain eq. (18). For the aforementioned case of pure perpendicular advection 
(i.e., 7 = 0), g(y) is a Gaussian function. A larger value of 7 implies a stronger power-law tail for y ^> 1. 

If we had allowed the boundary condition <?'(0) to be a nonzero constant, the resulting solutions for g(y) would 
have been asymmetric about fcy = 0. Naively, these solutions could be considered to be applicable to the relevant 
case of more power in outward propagating fluctuations than in inward propagating fluctuations. However, the solu- 
tions exhibit unphysically negative values of g throughout most of the "weaker" half of the distribution, and thus are 
probably inapplicable. 
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Fig. 1 . — Representative solutions for the dominant perpendicular cascade at r = 2 Rq : (a) outer-scale energy source 
function S(k±), normalized to its maximum value; (b) dissipation rate D(k±) = Do(k±R p ) 2 , where Do = 0.05 s _1 
(dotted line), 0.5 s _1 (dot-dashed line), and 5 s _1 (dashed line); (c) time-steady reduced power spectrum W±(k±) for 
the three dissipation rates above, and for zero dissipation (solid line). The proton gyroradius R p at this radius in the 
extended corona is assumed to be equal to 0.024 km. 
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Fig. 2. — Dispersion and energy density properties of kinetic Alfven waves as a function of k± : (a) real (solid line) 
and imaginary (dashed line) parts of the wave frequency, scaled by the MHD Alfven frequency, and the dimensionless 
energy parameter </> e (dotted line); (b) linearized energy density fractions — expressed as areas — from top to bottom: 
magnetic field, proton kinetic (i.e., velocity perturbation), electron kinetic, proton thermal (i.e., density perturbation), 
and electron thermal. All curves correspond to a constant value of fc|| Va/Q p = 10~ 3 . 
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Fig. 3. — Dispersion and energy density properties of nearly parallel-propagating Alfven waves (approaching ion 
cyclotron resonance) as a function of fcy. Curves in (a) and areas in (b) are similar to those in Figure 2, but the 
normalizations in (a) are different {see plot annotations). All curves correspond to a constant value of k±R p = 10~ 3 . 
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Fig. 4. — (a) Dimensionless damping rates \u)i tS /ui r \ for protons (unfilled contours) and electrons (filled contours) 
plotted one per decade between 9 x 10~ 5 and 9 x times the maximum values of \u)i iS /u r \ of 0.329 (protons) and 
0.231 (electrons). Both sets of contours go from low to high values with increasing wavenumber. Dotted lines show 
contours in 9, the angle between Bo and k. (b) Contours of VK(k), plotted one per 10 5 between 10 10 and 10 40 cm 5 s~ 2 
(solid lines) with the largest values at lower left. Critical balance curves also plotted for both the undamped (upper 
dashed line) and damped (lower dashed line) calculations of W^(k^). 
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Fig. 5. — Perpendicular proton heating rates, in units of a fiducial empirical rate Q emp (see text), plotted versus 
for a range of total wave amplitudes (SU /p) 1 ^ 2 (see labels on each curve). The hashed region denotes empirical rates 
from Li etal. (1999). 
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Fig. 6. — (a) Normalized proton heating rates plotted versus heliocentric radius, both for the turbulent cascade rate of 
Dmitruk et al. (2002) (dashed line) and for the derived ion/EPH diffusion coefficient (solid lines). From bottom to top, 
S = 2.5, 3, 3.5, 4. (b) Perpendicular proton temperatures derived with the ion/EPH heating rates plotted in (a) (solid 
lines), and with no imposed heating (i.e., adiabatic cooling only; dotted line). 
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Fig. 7. — Schematic representation of the major physical processes discussed in this paper. The relative amount 
of turbulent dissipation that directly heats protons and electrons depends, among other factors, on the ratio 0/<y. 
The nonlinear development of parallel electron beams into phase-space holes that can interact with protons is also 
illustrated. 



